Nonlocal pseudopotential energy density functional for orbital-free density functional theory

Orbital-free density functional theory (OF-DFT) is an electronic structure method with a low computational cost that scales linearly with the number of simulated atoms, making it suitable for large-scale material simulations. It is generally considered that OF-DFT strictly requires the use of local pseudopotentials, rather than orbital-dependent nonlocal pseudopotentials, for the calculation of electron-ion interaction energies, as no orbitals are available. This is unfortunate situation since the nonlocal pseudopotentials are known to give much better transferability and calculation accuracy than local ones. We report here the development of a theoretical scheme that allows the direct use of nonlocal pseudopotentials in OF-DFT. In this scheme, a nonlocal pseudopotential energy density functional is derived by the projection of nonlocal pseudopotential onto the non-interacting density matrix (instead of “orbitals”) that can be approximated explicitly as a functional of electron density. Our development defies the belief that nonlocal pseudopotentials are not applicable to OF-DFT, leading to the creation for an alternate theoretical framework of OF-DFT that works superior to the traditional approach.

A b initio calculations using Kohn-Sham (KS) density functional theory (DFT) 1,2 can accurately describe the fundamental properties of various materials. However, its computational cost scales with the cube of the number of electrons in the simulation cell, which poses a major challenge to large-scale simulations. In contrast, orbital-free (OF) DFT is inherent of the lower computational cost that scales linearly with the number of atoms in the system, as it relies only on the electron density and the use of KS orbitals is avoided. As a result, OF-DFT is successfully applied to large-scale simulations of systems with up to millions of atoms [3][4][5][6] .
The accuracy of OF-DFT simulations depends strongly on the quality of the non-interacting kinetic energy and the electron-ion (or electron-pseudocore) interaction energy employed in the simulations. Many approximate kinetic energy density functionals (KEDFs) have been proposed to evaluate the noninteracting kinetic energy in OF-DFT  . Their use in combination with local pseudopotentials [40][41][42][43][44] can achieve results that agree reasonably with those derived by KS-DFT, especially for main-group metals, III-V semiconductors 24,34,36,45 , and even systems with inhomogeneous electron density such as metal clusters and quantum dots 37,39 .
Unfortunately, the local pseudopotentials 28,31,32,40,41,43 used to evaluate the electron-ion interaction energy suffer from a lack of transferability 43 , as they fail to reproduce the correct scattering behavior of the all-electron potentials [46][47][48] . Overcoming the transferability problem requires a reliance on either all-electron potential or nonlocal pseudopotentials (NLPPs), which are widely used in orbital-based approaches. However, it is practically unfeasible to use the all-electron potential, as an accurate allelectron KEDF for OF-DFT calculations is not yet available 49,50 . Furthermore, the use of NLPPs 46,51 runs against conventional understanding, as no orbitals are available in the traditional framework of OF-DFT 41,42,44,52,53 .
A crucial nonlocal energy term with a set of angularmomentum-dependent energies has recently been added to OF-DFT calculations 54,55 in an effort to correct errors arising from the use of KEDFs and local pseudopotentials. This approach has successfully reproduced the bulk properties of several standard structures of Ti. However, special care must be taken when applying it to a wide range of practical simulations, as frozen onsite orbitals and empirically directed fitting parameters are part of the model 52 . There is substantial demand for a general approach to evaluate the electron-ion interaction energy using NLPPs in OF-DFT calculations. In this manuscript, we developed a theoretical scheme that allows the direct use of the NLPPs for the calculation of electron-ion interaction energy in OF-DFT, together with a specially designed theoretical framework of OF-DFT. This development leads to an OF-DFT calculation that gives a better transferability than the existing OF-DFT method based on local pseudopotentials.

Results and discussion
Nonlocal pseudopotential energy density functional. In general, the total energy density functional of OF-DFT can be expressed as: where ρ, T s , E H , E XC , E II , {R a }, and V loc are the electron density, KEDF, Hartree energy, exchange-correlation energy, ion-ion repulsion energy, the set of atomic positions, and local pseudopotential, respectively. To include the nonlocal electron-ion interactions, the total energy density functional of OF-DFT is reformulated as: where the total electron-ion interaction energy E EI [ρ] can be separated into two parts: a local part E loc [ρ] = ∫ρ(r)V loc (r)d 3 r and a nonlocal part E nl [ρ]. All of the terms in Eq.
The exact NLPP energy depends on the KS orbitals or the density matrix: where f i , V nl ðr 0 ; rÞ ¼ hr 0 jV nl jri, and γ s ðr; r 0 Þ ¼ ∑ i f i ψ i ðrÞψ Ã i ðr 0 Þ represent the occupation number of the ith KS orbital ψ i , the realspace representation of the nonlocal part pseudopotential, and the non-interacting density matrix, respectively. Considering that the density matrices γ s ½ρðr; r 0 Þ can be used to approximate the KEDFs 7,56,57 , an NLPP energy density functional (NLPPF) relying directly on the density matrix is proposed to evaluate the nonlocal electron-ion interaction energy. The nonlocal electron-ion interaction energy is then rewritten as a function of electron density where E a;lm KB ¼ ½ R ϕ aÃ lm ðrÞδV a l ðrÞϕ a lm ðrÞd 3 r À1 and χ a lm ðrÞ ¼ δV a l ðrÞ ϕ a lm ðrÞ. The terms ϕ a lm and δV a l are the atomic pseudo-wavefunction and the short-range pseudopotential corresponding to the lmth angular momentum of ath atom, respectively. γ s [ρ] denotes the density matrix as a function of electron distribution ρ. Although there is no exact analytic form available for the density matrix functional, a modified Gaussian (MG) 60 form derived from the second-order Taylor expansions of the density matrix 61 was employed to approximate the density matrix functional: where s ¼ jr À r 0 j and r ¼ ðr þ r 0 Þ=2. The second term in the square bracket is O(s 4 ) correction 60 , where A is an adjustable parameter. β(r) denotes the "local temperature" βðrÞ ¼ where t s (r) is the exact kinetic energy density defined as t s ðrÞ ¼ t KS s ðrÞ ∑ Occ: i¼1 1 8 j∇ρ i ðrÞj 2 =ρ i ðrÞ À 1 8 ∇ 2 ρðrÞ and ρ i (r) = |ψ i (r)| 2 is ith KS orbital's density (see Refs. 60,61 ). To remove the orbitaldependent problem in t s (r) and obtain a solely density-dependent form of the density matrix functional, the kinetic energy density is obtained directly from the integrand of KEDFs to replace the exact one: t s (r) ≈ t s [ρ](r). The widely used Wang-Teter (WT) KEDF 28 is chosen as an exemplary case, and t s [ρ](r) can be expressed as t WT s ½ρðrÞ ¼ where ω WT ðr; r 0 Þ is the kernel of WT functional. The Supplementary Notes give the details of the kinetic energy densities obtained from WT and Xu-Wang-Ma 38 KEDFs. The direct numerical evaluations of ρð rÞ and βð rÞ at the average position r are very complicated. They are therefore approximated using q-mean "nonlocal density" ρ q ðr; r 0 Þ ¼ ρ q ðrÞþρ q ðr 0 Þ 2 h i 1=q and two-point average temperature βðr; r 0 Þ ¼ ½βðrÞ þ βðr 0 Þ=2 for systems with slowly varying electron densities. The density matrix functional of Eq. (6) can then be reformulated as Combining Eqs. (4), (5), and (8) gives the NLPPF as where the integral domain Ω a is the ath ionic core region. Owing to the short-range nature of fχ a lm ðrÞg, the computational cost of Eq. (9) scales linearly O½cN a with the number of atoms (N a ), where c can be regarded as a constant derived from the double integral within the near-core region. Within the NLPPF scheme, a new theoretical framework of OF-DFT has been built and implemented in ATLAS 5,64 . The further computational details are provided in Methods Section. The parameters of pseudopotential and NLPPF are presented in Supplementary Tables 1 and 2, respectively.
Computational accuracy of NLPPF scheme. We first applied this scheme for OF-DFT calculations of Li, Mg, and Cs within hexagonal-close-packed (HCP), face-centered cubic (FCC), bodycentered cubic (BCC), simple cubic (SC), and cubic diamond structures. For each structure, 13 energy-volume points were calculated by expanding and compressing the approximate equilibrium volume by up to 20%, and the bulk properties (the equilibrium cell volume V 0 , bulk modulus B 0 , and the relative energy E R with respect to HCP structure) were determined by fitting the energy-volume curve against Murnaghan's equation of state 65 . The comparison of the results obtained by OF-DFT using both local pseudopotentials (BLPS 41 and OEPP 43 ) and NLPPs against those calculated by KS-DFT using the projector augmented-wave (PAW) 66 method are presented in Supplementary Table 3. For Li and Mg solids, the OF-DFT calculations within both local pseudopotentials and NLPPs give reasonable predictions of V 0 , B 0 , and E R , which are comparable with the KS-DFT results. It is noteworthy that our scheme shows an improvement over the local pseudopotentials of OEPP for bulk Cs. The accurate bulk properties of Li/Mg/Cs obtained by the current scheme demonstrate its valid applicability to simple metallic solids.
Further assessment of the accuracy of our scheme was demonstrated by molecular dynamics calculations for Li-Mg alloy. The calculations used a canonical ensemble (at 1000 K) in a supercell containing 108 atoms (Li 54 Mg 54 ). The calculated pair distribution functions g(r) for Li-Mg alloy are shown in Fig. 1. Overall, the predicted shapes and peak neighbors of pair distribution functions by OF-DFT within NLPPF match the  (Fig. 1b-d) calculated by OF-DFT within the NLPPF being almost identical to the KS-DFT calculations, which are superior to those obtained by the local pseudopotentials.
Transferability of NLPPF scheme. To demonstrate the transferability of our scheme, we randomly generated 50 structures of Li systems using CALYPSO 67,68 . The total energies of these structures were calculated by OF-DFT and KS-DFT. The comparisons of total energy relative to the HCP structure (E R ) are shown in Fig. 2. The orderings of energy are well captured by OF-DFT within NLPPF. The relative energies of different phases are overall well reproduced and in reasonable agreement with the KS-DFT results. The least-square fitting lines of WT-NLPPF are generally closer to the KS-PAW results than those from the local pseudopotentials. For example, the mean error of E R for Li systems obtained by OF-DFT within the WT-NLPPF is 45 meV/ atom, which is lower than that within BLPS (73 meV/atom), or OEPP (201 meV/atom). Therefore, this framework of OF-DFT with improved transferability is superior to the traditional one.
Previous studies have shown that OF-DFT with local pseudopotentials can be applied to most s-and p-block metals. However, OF-DFT simulation using the local pseudopotential OEPP shows unacceptable errors for various crystalline phases of Be (Fig. 3): the curves of energy with respect to volume for HCP and FCC structures show the total energy monotonically decreasing with increasing volume. In contrast, the curves with clear minima predicted by OF-DFT within NLPPF agree well with those produced by KS-DFT. These findings indicate the significant superiority of our proposed framework over the conventional one.
Note that the bulk properties (e.g., equilibrium volume, bulk modulus, and relative energy) calculated by OF-DFT within NLPPF reproduce the results of KS-DFT for Li and Mg almost exactly (Table 1), considering the maximal deviation of E R is within 32 meV/atom. However, there are some discrepancies for crystalline phases of Be: in particular, the deviation of E R for the SC structure is larger than 400 meV/atom. To explore the causes of these discrepancies, we estimated the errors of the kinetic energy density of the WT-KEDF with respect to the KS kinetic energy density along the [100] and [111] directions in the SC structures of Li, Mg, and Be (Fig. 4). The kinetic energy density of KS-DFT is clearly reproduced accurately by the WT-KEDF for Li and Mg with slow variations of electron densities. However, it is seriously underestimated for Be, in which the electron distribution rapidly varies in the near-core region. Therefore, we believe that errors in the kinetic energy densities for Be lead to the discrepancy in its bulk properties obtained by the framework of   OF-DFT within the NLPPF. The findings are fairly consistent with our expectation that the performance of the NLPPF relies strongly on the accuracy of kinetic energy density, as manifested by Eqs. (7) and (8).
Due to the existence of significant differences in kinetic energy density between WT and the exact one for Cd systems including localized d-channel electrons (see Supplementary Fig. 1), the NLPPF using WT cannot be applicable to investigate Cd systems with d-channel electrons. Therefore, the NLPP of Cd was constructed without d-channel electrons for the additional calculations. As listed in Table 2, the bulk properties predicted by the OF-DFT within NLPPF framework agree fairly well with the predictions by KS-DFT using the same NLPP, whereas WT-OEPP gives serious discrepancies compared with KS-NLPP results in all bulk properties. Although the calculations using NLPP without d-channel electrons cannot give the accurate bulk properties for Cd systems (Table 2), OF-DFT within NLPPF works superior to that within OEPP. Overall, it can be expected this NLPPF scheme using accurate KEDF and its kinetic energy density can be applied to the systems including localized electrons, such as the transition metals or covalent systems.
Computational efficiency of NLPPF scheme. To assess the computational efficiency of the current scheme, static simulations of Cs BCC supercells containing 128-16,000 atoms were performed by OF-DFT within NLPPF. The total wall time of the single point energy calculations is plotted with respect to the number of atoms in Supplementary Fig. 2. The computational cost of this framework clearly scales linearly with the number of atoms in the simulation cell, in sharp contrast to the cubic scaling of KS-DFT. This shows that the OF-DFT within NLPPF is potentially applicable to the simulation of large-scale systems containing millions of atoms.
In summary, we proposed an NLPPF scheme that allows the direct use of NLPPs in OF-DFT calculations. The static and dynamic properties of s-and p-block metals calculated within this scheme agree well with KS-DFT predictions and show significant improvements in the computational accuracy and transferability over conventional OF-DFT with local pseudopotentials. With this work, we defy the conventional wisdom of orbital-dependent NLPPs being incompatible with OF-DFT, leading to the creation of an alternative framework of OF-DFT, which opens up new avenues for further development of the theory.

Methods
Pseudopotential generations. The Troullier-Martins NLPPs 59 are generated by the FHI98PP 47 code for all considered systems [see Supplementary Table 1] and the p-channel of the NLPPs is used as the local pseudopotential of V loc (r) in OF-DFT.
Numerical calculations. The KS-DFT calculations using the PAW 66 and NLPP are performed by VASP 69,70 and ARES packages 71 , respectively. The k-point meshes are generated using the Monkhorst-Pack method 72 with the k-spacing of 0.10 Å −1 . The kinetic energy cutoff is 500 eV for all the simulations using VASP. The OF-DFT calculations are carried out by ATLAS 5,64 using WT 28 as KEDF, and the corresponding kinetic energy density is used to construct the NLPPF. The generalized gradient approximation with the form of Perdew-Burke-Ernzerhof 73  Molecular dynamics. The molecular dynamic simulations of Li-Mg alloy are performed in the canonical ensemble (at 1000 K) applying the Nosé-Hoover thermostat 74,75 simulations up to 10 ps (0.5 fs/step), with the first 10,000 steps for equilibrating the system. The data for further analysis were collected from the subsequent 10,000 steps.

Data availability
The authors declare that the main data supporting the findings of this study are contained within the paper and its associated Supplementary Information. All other relevant data are available from the corresponding authors upon reasonable request.